Matrix elements relevant for AI = 1/2 rule and e' / e from 
: Lattice QCD with staggered fermions 

On 



c3 



O 



00 
» 

Q 



D. Pekurovsky and G. Kilcup 



Department of Physics, 
the Ohio State University, 
<N : 174 W. 18th Ave., Columbus OH 43210, USA 

On 



Abstract 



We perform a study of matrix elements relevant for the AI = 1/2 rule and the direct 
Q\ • CP-violation parameter e'/e from first principles by computer simulation in Lattice QCD. 

We use staggered (Kogut-Susskind) fermions, and employ the chiral perturbation theory 
method for studying K° — ► irir decays. Having obtained a reasonable statistical accuracy, 
we observe an enhancement of the AI = 1/2 amplitude, consistent with experiment 
within our large systematic errors. Finite volume and quenching effects have been studied 
and were found small compared to noise. The estimates of e'/e are hindered by large 
uncertainties associated with operator matching. In this paper we explain the simulation 
^ ■ method, present the results and address the systematic uncertainties. 
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1 Introduction 



In those areas of particle phenomenology which require addressing non-perturbative effects, 
Lattice QCD plays an increasingly significant role, being a first-principles method. The rapid 
advances in computational performance as well as algorithmic techniques are allowing for better 
control of various errors associated with lattice calculations. 

In this paper we address the phenomenology of K° — > tiit decays. One of the long-standing 
puzzles is the "A/ =1/2 rule" , which is the observation that the transition channel with isospin 
changing by 1/2 is enhanced 22 times with respect to transitions with isospin changing by 3/2. 
The strong interactions are essential for explaining this effect within the Standard Model. 
Since the energy scales involved in these decays are rather small, computations in quantum 
chromodynamics (QCD) have to be done using a non-perturbative method such as Lattice 
QCD. Namely, Lattice QCD is used to calculate hadronic matrix elements of the operators 
appearing in the effective weak Hamiltonian. 

There have been so far several other attempts to study matrix elements of the operators 
relevant for AI = 1/2 rule on the lattice J| [|, |6|], but they fell short of desired accuracy. In 
addition, several groups |7|, §] have studied matrix elements (7r + 7r°|Oj|i ; r + ), which describe only 
AI = 3/2, not AI = 1/2 transition. In the present simulation, the statistics is finally under 
control for AI = 1/2 amplitude. 

Our main work is in calculating matrix elements (ix + \Oi\K + ) and (0\Oi\K°) for all basis 
operators (introduced in Sec. [2.ip . This is enough to recover matrix elements (7nr\Oi\K°) using 
chiral perturbation theory in the lowest order, although this procedure suffers from uncertainties 
arising from ignoring higher orders (in particular, final state interactions). The latter matrix 
elements are an essential part of the phenomenological expressions for AI = 1/2 and AI = 3/2 
amplitudes, as well as e'/e. The ratio of the amplitudes computed in this way confirms signifi- 
cant enhancement of AI = 1/2 channel, although systematic uncertainties preclude a definite 
answer. 

In addition, we address a related issue of e'/e - the direct CP- violation parameter in the 
neutral kaon system. As of the day of writing, the experimental data are somewhat ambiguous 
about this parameter: the group at CERN (NA48) reports Re(e'/e) = (23 ± 7) x 10" 4 , 
while the Fermilab group (E731) M has found ~Re(e'/e) = (7.4 ± 6.0) x 10~ 4 . There is a hope 
that the discrepancy between the two reports will soon be removed in a new generation of 
experiments. 

On the theoretical side, the progress in estimating e'/e in the Standard Model is largely 
slowed down by the unknown matrix elements [12] of the appropriate operators. The previous 
attempts |3|, f|, || to compute them on the lattice did not take into account operator matching. 
In this work we repeat this calculation with better statistics and better investigation of sys- 
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tematic uncertainties. We are using perturbative operator matching. In some cases it does not 
work, so we explore alternatives and come up with a partially non-perturbative renormaliza- 
tion procedure. The associated errors are estimated to be large. This is currently the biggest 
stumbling block in computing e'/s. 

The paper is structured as follows. In the Section || we show the context of our calculations, 
define the quantities we are looking after and discuss a number of theoretical points relevant for 
the calculation. Section § discusses issues pertaining to the lattice simulation. In Section |] we 
present the results and discuss systematic errors for AI = 1/2 rule amplitudes. In Section | we 
explain how the operator matching problem together with other systematic errors preclude a 
reliable calculation of e'/s, and give our best estimates for this quantity in Section Section [7] 
contains the conclusion. In the Appendix we give details about the quark operators and sources, 
and provide explicit expressions for all contractions and matrix elements for reference purposes. 



2 Theoretical framework 



2.1 Framework and definitions 



The standard approach to describe the problems in question is to use the Operator Product 
Expansion at the My/ scale and use the Renormalization Group equations to translate the 
effective weak theory to more convenient scales (/i ~ 2-4 GeV). At these scales the effective 



Hamiltonian for K — > ixix decays is the following linear superposition |T2 
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H w = -tJk^EN)") + TVi(ri\OiQi) , (i) 



where Zi and yi are Wilson coefficients (currently known at two-loop order), r = —V t dV t * s /V u dV*. 
and Oi are basis of four-fermions operators defined as follows: 
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(11) 

Here a and (3 are color indices, e q is quark electric charge, and summation is done over all light 
quarks. 

Isospin amplitudes are defined as 

A , 2 e iSo ' 2 EE ((nn) I= o, 2 \H w \K }, (12) 
where £0,2 are the final state interaction phases of the two channels. Experimentally 

u = ReA /ReA 2 ~ 22 . (13) 
Direct CP violation parameter e' is defined in terms of imaginary parts of these amplitudes: 

, Iim4n — tdlmAo u /o,a n / \ 

V2uReA 

Experiments are measuring the quantity Ree'/e, which is given by 

Re 7^27^ ImA ' |n °- n J- (15) 

where 

n = '£y i ((™)i =0 \o? >) \K )(i-n rtW ) (16) 

i 

n 2 = ^((TTTr^lOfV) (17) 

i 

with ImA t ee ImV^VJ*, and where f^+r,' ~ 0.25 ± 0.05 takes into account the effect of isospin 
breaking in quark masses (m u 7^ m^). Oj- ^ and o| 2 ^ are isospin and 2 parts of the basis 
operators. Their expressions are given in the Appendix for completeness. 
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2.2 Treatment of charm quark 



The effective Hamiltonian given above is obtained in the continuum theory in which the top, 
bottom and charm quarks are integrated out. (In particular, the summation in Eqs. (fiHTTD 
is done over u, d and s quarks.) This makes sense only when the scale /i is sufficiently low 
compared to the charm quark mass. As mentioned in Ref. |§, at scales comparable to m c higher- 
dimensional operators can contribute considerably. Then one should consider an expanded set 
of operators including those containing the charm quark. Lattice treatment of the charm quark 
is possible but in practice quite limited, for example by having to work at much smaller lattice 
spacings and having a more complicated set of operators and contractions. Therefore we have 
opted to work in the effective theory in which the charm quark is integrated out. Since we 
typically use /i ~ 2 GeV in our simulations, this falls into a dangerous region. We hope that 
the effects of higher-dimensional operators can still be neglected, but strictly speaking this issue 
should be separately investigated. 



2.3 Calculating (7nr\Oi\K 



As was shown by Martinelli and Testa |TJ]], two-particle hadronic states are very difficult to 
construct on the lattice (and in general, in any Euclidean description). We have to use an alter- 
native procedure to calculate the matrix elements appearing in Eqs. (O'lC'O)- We choose the 
method (Tl[| in which the lowest-order chiral perturbation theory is used to relate (iT7T\Oi\K ) 
to matrix elements involving one-particle states: 



(7T 7T 



\Oi\K°) 



m K — mi 



f 



-1 



(n + \O l \K + ) 
(0\Oi\K°) 



m s + m d 
[P-k 'Pk)1 : o 



f 



m d )5, 



(19) 
(20) 



where / is the lowest-order pseudoscalar decay constant. The masses in the first of these 
formulae are the physical meson masses, while the quark masses and the momenta in the 
second and third formulae are meant to be from actual simulations on the lattice (done with 
unphysical masses). These relationships ignore higher order terms in the chiral expansion, 
most importantly the final state interactions. Therefore this method suffers from a significant 
uncertainty. Golterman and Leung |H| have computed one-loop correction for AI = 3/2 



amplitude in chiral perturbation theory. They find this correction can be large, up to 30% or 
60%, depending on the values of unknown contact terms and the cut-off. 
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3 Lattice techniques 



3.1 Mixing with lower-dimensional operators. 

Eqs. ([18|-|2"0|) handle unphysical s <-> d mixing in (7r + |Oj|/i + ) by subtracting the unphysical 
part proportional to (0\Oi\K°). This is equivalent to subtracting the operator 

O sub = (m d + m s )sd + (m d - m s )sj 5 d . (21) 

As shown by Kilcup, Sharpe et al. in Refs. || [K|, these statements are also true on the lattice 



if one uses staggered fermions. A number of Ward identities discussed in these references 
show that lattice formulation with staggered fermions retains the essential chiral properties 
of the continuum theory. In particular, O su b defined in Eq. [H] is the only lower-dimensional 
operator appears in mixing with the basis operators. (Lower-dimensional operators have to 
be subtracted non-perturbatively since they are multiplied by powers of a -1 .) We employ the 
non-perturbative procedure suggested in Ref. (TT 



(n + n-\O t \K°) = (tt+IO, - a t O suh \K + ) ■ ™" K ™* , (22) 

(Pn-PK)f 

where ati are found from 

= (0\O l -a t O sub \K°). (23) 

This procedure is equivalent to the lattice version of Eqs. (pT8|-pO|) and allows subtraction times- 
lice by timeslice. 

Throughout our simulation we use only degenerate mesons, i.e. m s = = m u . Since only 
negative parity part of O su b contributes in Eq. (|23D , one naively expects infinity when calculating 
OLi. However, matrix elements (0\Oi\K°) of all basis operators vanish when m s = rrid due to 
invariance of both the Lagrangian and all the operators in question under the CPS symmetry, 
which is defined as the CP symmetry combined with interchange of s and d quarks. Thus 
calculation of requires taking the first derivative of (0\Oi\K°) with respect to (m^ — m s ). In 
order to evaluate the first derivative numerically, we insert another fermion matrix inversion in 
turn into all propagators involving the strange quark. Detailed expressions for all contractions 
are given in the Appendix. 

3.2 Diagrams to be computed 



According to Eqs. (|22| , p3|) , we need to compute three diagrams involving four- fermion operators 
(shown in Fig. [I]) and a couple of bilinear contractions. The "eight" contraction type (Fig. [L]a) 
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(a) 



(b) 




Figure 1: Five diagrams types needed to be computed: (a) "Eight"; (b) "Eye"; (c) "Annihila- 
tion"; (d) "Subtraction"; (e) two-point function. 



6 



is relatively cheap to compute. It is the only contraction needed for the AI = 3/2 amplitude. 
The "eye" and "annihilation" diagrams (Fig. [l|b and |l|c) are much more expensive since they 
involve calculation of propagators from every point in space-time. 



3.3 Lattice parameters and other details 

The parameters of simulation are listed in the Table [l]. We use periodic boundary conditions 
in both space and time. Our main "reference" ensemble is a set of quenched configurations at 
(3 = Q/g 2 = 6.0 (Qi). In addition, we use an ensemble with a larger lattice volume (Q2), an 
ensemble with (3 = 6.2 (Q 3 ) for checking the lattice spacing dependence, and an ensemble with 
2 dynamical flavors (m = 0.01) generated by the Columbia group, used for checking the impact 
of quenching. The ensembles were obtained using 4 sweeps of SU(2) overrelaxed and 1 sweep 
of 577(2) heatbath algorithm^. The configurations were separated by 1000 sweeps, where one 
sweep includes three SU(2) subgroups updates. 



Table 1: Simulation parameters 



Ensemble 


N f 


P 




Size 






L, fm 


Number of 


Quark masses 


name 
















configurations 


used 


Qi 





6.0 


16 3 


x (32 


x 


4) 


1.6 


216 


0.01 — 0.05 


Q2 





6.0 


32 3 


x (64 


x 


2) 


3.2 


26 


0.01 — 0.05 


Qs 





6.2 


24 3 


x (48 


x 


4) 


1.7 


26 


0.005 — 0.03 


D 


2 


5.7 


16 3 


x (32 


X 


4) 


1.6 


83 


0.01 — 0.05 



We use the standard staggered fermion action. Fermion matrices are inverted by conjugate 
gradient. Jackknife is used for statistical analysis. 

As explained below, we have extended the lattice 4 timesf] in time dimension by copying 
gauge links. This is done in order to get rid of excited states contamination and wrap-around 
effects. 

The lattice spacing values for quenched ensembles were obtained by performing a fit in the 
form of the asymptotic scaling to the quenched data of p meson mass given elsewhere |TT |. 



Lattice spacing for the dynamical ensemble is also set by the p mass fl6 



Some other technicalities are as follows. We work in the two flavor formalism. We use local 
wall sources that create pseudoscalar mesons at rest. (Smearing did not have a substantial 
effect.) The mesons are degenerate (m s = rrid = m u , m n = tuk). We use staggered fermions and 



1 except for the dynamical set which was obtained by R- algorithm ]iq| 
2 for all ensemble except the biggest volume, which we extend 2 times. 
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work with gauge-invariant operators, since the gauge symmetry enables significant reduction 
of the list of possible mixing operators. The staggered flavour structure is assigned depending 
on the contraction type. Our operators are tadpole-improved. This serves to 'improve" the 
perturbative expansion at a later stage when we match the lattice and continuum operators. 
For calculating fermion loops we employ the U(l) pseudofermion stochastic estimator. More 
details and explanation of some of these terms can be found in the Appendix. 




Figure 2: The general setup of calculation of (7r + \Oi\K + ) (without loss of generality, an "eight" 
contraction is shown). The kaon source is at the timeslice 0, while the pion sink is at the 
timeslice T. The operator is inserted at a variable time t. The result of this contraction is 
proportional to the product of two exponentials shown in the figure. 



3.4 Setup for calculating matrix elements of four-fermion 
operators 

Consider the setup for calculation of (ji + \Oi\K + ) . Kaons are created at t = 0, the operators 
are inserted at a variable time t, and the pion sink is located at the time T (see Fig. ^), where 
T is sufficiently large. In principle, a number of states with pseudoscalar quantum numbers can 
be created by the kaon source. Each state's contribution is proportional to \/~Ze~ ra \ t \ so the 
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lightest state (kaon) dominates at large enough t. Analogously, states annihilated by the sink 
contribute proportionally to \[Ze^ m}[T ~ t \ which is dominated by the pion. 

In this work kaon and pion have equal mass. In the middle of the lattice, where t is far 
enough from both and T, we expect to see a plateau, corresponding to Ze~ m7rT (tt\0\K). This 
plateau is our working region (see Fig. |]). 

As concerns the kaon annihilation matrix elements (0\Oi\K°), we only need their ratio to 
(Ols75dl.fr ), in which the factors \/~Ze~ mt cancel. Indeed, we observe a rather steady plateau 
(Fig. |). 




Figure 3: B ratio is formed by dividing the four-fermion matrix element by the product of two- 
point functions, typically involving or P bilinears. All the operators involved are inserted at 
the same timeslice t, and summation is done over spatial volume. The external meson sources 
are located at timeslices and T, both in the numerator and the denominator. This enables 
cancellation of some common factors. 

3.5 B ratios 

It has become conventional to express the results for matrix elements in terms of so-called B 
ratios, which are the ratios of desired four-fermion matrix elements to their values obtained 
by vacuum saturation approximation (VSA). For example, the B ratios of operators O2 and 
O4 are formed by dividing the full matrix element by the product of axial-current two-point 
functions (Fig. |3|). We expect the denominator to form a plateau in the middle of the lattice, 
equal to Ze~ m ^ T (7t|t4 m |0) ■ (Ol-A^I-fQ, where A^ are the axial vector currents with appropriate 
flavor quantum numbers for kaon and pion. The factor Ze~ m7lT cancels, leaving the desirable 
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ratio (ir\0\K) / ((^1^10) • (0|A M |iT)). Apart from common normalization factors, a number 
of systematic uncertainties also tend to cancel in this ratio, including the uncertainty in the 
lattice spacing, quenching and in some cases perturbative correction uncertainty. Therefore, it 
is sometimes reasonable to give lattice answers in terms of the B ratios. 

However, eventually the physical matrix element needs to be reconstructed by using the 
known experimental parameters (namely fx) to compute VSA. In some cases, such as for 
operators 5 — 8 , the VSA itself is known very imprecisely due to the failure of perturbative 
matching (see Sec. ||). Then it is more reasonable to give answers in terms of matrix elements 
in physical units. We have adopted the strategy of expressing all matrix elements in units of 
(7r|A M |0) (0|A^|iT) = (fK tt ) 2 m 2 M at an intermediate stage, and using pre-computed f^ tt at the 
given meson mass to convert to physical units. This method is sensitive to the choice of the 
lattice spacing. 




10 20 30 40 50 

t 



Figure 4: An example of the signal we get for one of the B ratios (in this case, for the "eye" 
part of the O2 operator on Qi ensemble). The wall sources are at t — 1 and t = 49. We see that 
the excited states quickly disappear and a stable, well-distinguished plateau is observed. We 
perform jackknife averaging in the range of t from 12 to 37 (shown with the horizontal lines). 
It is important to confirm the existence of the plateau for reliability of the results. 

It is very important to check that the time distance between the kaon and pion sources T 
is large enough so that the excited states do not contribute. That is, the plateau in the middle 
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Figure 5: An example of the signal for (O|02|J£' ) / l( m d ~ rn s ) {^\sj5d\K )] on Qi ensemble. 
The kaon source is at t — 1. We average over the range of t from 5 to 12 (shown with horizontal 
lines) . 

of the lattice should be sufficiently flat, and the B ratios should not depend on T. We have 
found that in order to satisfy this requirement the lattice has to be artificially extended in time 
direction by using a number of copies of the gauge links (4 in the case of the small volume 
lattices, 2 otherwise). We are using T = 72 for Q 3 {(3 = 6.2) ensemble, and T = 48 for the rest. 
An example of a plateau that we obtain with this choice of T is shown in Fig. ^. To read off 
the result, we average over the whole extension of the plateau, and use jackknife to estimate 
the statistical error in this average. 

4 AI = 1/2 rule results 

Using the data obtained for matrix elements of basis operators, in this section we report numer- 
ical results for ReA and ReA 2 amplitudes as well as their ratio. We discuss these amplitudes 
separately since the statistics for KeA 2 is much better and the continuum limit extrapolation 
is possible. 
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4.1 Re^ results 

The expression for KeA 2 can be written as 

ReA 2 = ^ V ud V* s z+0m){O 2 )2, (24) 

where z + (fi) is a Wilson coefficient and 

(0 2 ) 2 = ((nn) I=2 \0 { 2 2) \K). (25) 

Here 

{2) = {2) = l[(^(l-75V)(«7 M (l-75)rf) 

+(s 7m (1 - 7B )d)(M7 M (l - 7s)w) - (*7„(1 - 75)^(^(1 - 7s)d)- (26) 
In the lowest order chiral perturbation theory the matrix element (02)2 can be expressed as 

ml -ml ( 7 r+\O i 2 2) \K+) , x 

(0 2 ) 2 = V2 K - n V 1 2 ' ^. 27 

/ 77i 

The latter matrix element involves only "eight" diagrams. Moreover, in the limit of preserved 
5 , f/(3)fl av0 r symmetry it is directly related |20| to parameter Bk (which is the B ratio of the 
neutral kaon mixing operator Ok = {s^id) (stlG?)), so that 



in - 4 v2 ■»//,--///- 2 P ,,,, 
(°^--g ^p-Aatt 5 - ^ 



2 2 
mi — m: 



Parameter Bk is rather well studied (for example, by Kilcup, Pekurovsky |ZTJ and JLQCD 



collaboration f22|). Quenched chiral perturbation theory [|6|] predicts the chiral behaviour of 



the form Bk = a + bm\ + c m\ logm K , which fits the data well (see Fig. ^) and yields a finite 
non-zero value in the chiral limit. Note that Rev4 2 is proportional to the combination BkJ^ 
and since both multipliers have a significant dependence on the meson mass (Figs. || and 0), 
KeA 2 is very sensitive to that mass. Fig. |8] shows KeA 2 data for the dynamical ensemble, 
based on B K values we have reported elsewhere pi| . Which meson mass should be used 
to read off the result becomes an open question. If known, the higher order chiral terms 
would remove this ambiguity. Forced to make a choice, we extrapolate to M 2 = (m K + m 2 r )/2. 
Using our data for Bk in quenched QCD and taking the continuum limit we obtain: KeA 2 = 
(1.7 ±0.1) ■ 10 -8 GeV, where the error is only statistical, to be compared with the experimental 
result ReA 2 = 1.23 • 1(T 8 GeV. 
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Figure 6: Parameter Bx in NDR MS scheme at 2 GeV on the dynamical ensemble vs. the 
meson mass squared. The fit is of the form Bk = a + bm 2 K + cm 2 K log m 2 K . The vertical line 
here and in the other plots below marks the physical kaon mass. 




Figure 7: Pseudoscalar decay constant (F n = 93 MeV experimentally) on the dynamical en- 
semble vs. meson mass squared. The fit is of the same form as Bk- 
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Figure 8: Re^ for the dynamical ensemble. The fit is of the same form as Bk- The horizontal 
line is the experimental value of 1.23 GeV. 

Higher order chiral terms (including the meson mass dependence) are the largest systematic 
error in this determination. According to Golterman and Leung ]Tj|, one- loop corrections 
in (quenched) chiral perturbation theory are expected to be as large as 30% or 60%. Other 
uncertainties (from lattice spacing determination, from perturbative operator matching and 
from using finite lattice volume) are much smaller. 

4.2 Re^4 results 

Using Eqs. (|2"2] , |2"3~D , ReA can be expressed asQ 

ReA = ^V ud V: s ^^ £ ztRt, (29) 

3 In our normalization Re^lo = 27.2 • 10 _s . 
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where Zi are Wilson coefficients and 



tii — 5 • 

mr 

The subscript V indicates that these matrix elements already include subtraction of 
(7i + \O su b\K + ). All contraction types are needed, including the expensive "eyes" and "anni- 
hilations". Of 1 are isospin parts of operators Oi (given in the Appendix for completeness). 
For example, 

{ ? = ^(^(l-75)rf)^7 /1 (l-75)w)-^(^(l-75)w)(^(l-7 5 )rf) 

+ 1(5^(1 - l5 )d)(dY(l-l 5 )d) (30) 

Of = 1(81,(1- l5 )u)(uY(l- ls)d) - 1(^(1 - 76 )d)(!T7' i (l-75)M) 

+ i(^(l-75)rf)(^(l-75)rf) (31) 

The results for quenched /J = 6.0 and (3 = 6.2 ensembles are shown in Fig. || Dependence 
on the meson mass is small, so there is no big ambiguity about the mass prescription as in the 



ReA 2 case. Some lattice spacing dependence may be present (Fig. |10|), although the statistics 
for (3 = 6.2 ensemble is too low at this moment. 

The effect of the final state interactions (contained in the higher order terms of the chiral 
perturbation theory) is likely to be large. This is the biggest and most poorly estimated 
uncertainty. 

An operator matching uncertainty arises due to mixing of 2 with 0$ operator through 



penguin diagrams in lattice perturbation theory. This is explained in the Section |5.1] . We 
estimate this uncertainty at 20% for all ensembles. 

As for other uncertainties, we have checked the lattice volume dependence by comparing 
ensembles Qi and Q 2 (1.6 and 3.2 fm at (3 — 6.0). The dependence was found to be small, so 
we consider (1.6 fm) 3 as a volume large enough to hold the system. We have also checked the 
effect of quenching and found it to be small compared to noise (see Fig. |TT]) . 



The breakdown of contributions of various basis operators to ReAo is shown in Fig. By 
far, O2 plays the most important role, whereas penguins have only a small influence. 
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Figure 9: ReA for quenched ensembles plotted against the meson mass squared. The upper 
group of points is for ensembles Qi and Q 2 , while the lower group is for Q 3 . Only statistical 
errors are shown. 
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Figure 10: KeA for quenched ensembles plotted against lattice spacing squared. The horizontal 
line shows the experimental result of 27.2 • 1CT 8 GeV. Only statistical errors are shown. 
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Figure 11: Comparison of quenched (Qi) and dynamical results for ReA at comparable lattice 
spacings. 
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Figure 12: Contribution of various operators to ReA . 



4.3 Amplitude ratio 

Shown in Fig. [13] is the ratio ReA 2 /ReA as directly computed on the lattice for quenched and 
dynamical data sets. The data exhibit strong dependence on the meson mass, due to Re^ 
chiral behaviour (compare with Fig. |8|). 

Within our errors the results seem to confirm, indeed, the common belief that most of 
the AI = 1/2 enhancement comes from the "eye" and "annihilation" diagrams. The exact 
amount of enhancement is broadly consistent with experiment while being subject to consid- 
erable uncertainty due to higher-order chiral terms. Other systematic errors are the same as 
those described in the previous Subsection. 

5 Operator matching 

As mentioned before, we have computed the matrix elements of all relevant operators with an 
acceptable statistical accuracy. These are regulated in the lattice renormalization scheme. To 
get physical results, operators need to be matched to the same scheme in which the Wilson 
coefficients were computed in the continuum, namely MS NDR. While perturbative matching 
works quite well for ReAo and ReA2, it seems to break down severely for matching operators 
relevant for e'/e. 
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Figure 13: The ratio Re^/Re^o versus the meson mass squared for quenched and dynamical 
ensembles. Ensembles Q\ and D have comparable lattice spacings. The dynamical ensemble 
data were used for the fit. The big slope of the fit line is accounted for by the mass dependence 
of Rev4 2 . The horizontal line shows the experimental value of 1/22. The error bars show only 
the statistical errors obtained by jackknife. 



5.1 Perturbative matching and Re^o 

Conventionally, lattice and continuum operators are matched using lattice perturbation theory: 
Or\q*) = O l t l + ln(,f a) + C^O? + 0(g 4 ) + 0(a n ), (32) 

where 7^ is the one-loop anomalous dimension matrix (the same in the continuum and on the 
lattice), and Cjj are finite coefficients calculated in one- loop lattice perturbation theory |25| 24 1. 
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We use the "horizontal matching" procedure [[HJ, whereby the same coupling constant as in the 
continuum (gjjg) is used. The operators are matched at an intermediate scale q* and evolved 
using the continuum renormalization group equations to the reference scale /i, which we take 
to be 2 GeV. 

In calculation of R,eA and KeA 2 , the main contributions come from left-left operators. One- 
loop renormalization factors for such (gauge-invariant) operators were computed by Ishizuka 
and Shizawa |25| (for current- current diagrams) and by Patel and Sharpe [24[ (for penguins). 



These factors are fairly small, so at the first glance the perturbation theory seems to work well, 
in contrast to the case of left-right operators essential for estimating e'/e, as described below. 
However, even in the case of KeA there is a certain ambiguity due to mixing of 2 operator with 
Oq through penguin diagrams. The matrix element of Oq is rather large, so it heavily influences 
(O2) in spite of the small mixing coefficient. Operator 0§ receives enormous renormalization 
corrections in the first order, as discussed below. Therefore, there is an ambiguity as to whether 
the mixing should be evaluated with renormalized or bare 0§. Equivalently, the higher-order 
diagrams (such as Fig. |T6| b and |T6|d) may be quite important. 



In order to estimate the uncertainty of neglecting higher-order diagrams, we evaluate the 
mixing with 0$ renormalized by the partially non-perturbative procedure described below, and 
compare with results obtained by evaluating mixing with bare Oq. The first method amounts 
to resummation of those higher-order diagrams belonging to type (b) in Fig. ITB, while the 



second method ignores all higher-than-one-order corrections. Results quoted in the previous 
Section were obtained by the first method, which is also close to using tree-level matching. The 
second method would produce values of ReA lower by about 20%. Thus we consider 20% a 
conservative estimate of the matching uncertainty. 

In calculating e'/e the operator matching issue becomes a much more serious obstacle as 
explained below. 

5.2 Problems with perturbative matching 

The value of e'/e depends on a number of subtle cancellations between matrix elements. In 
particular, Oq and 0§ have been so far considered the most important operators whose con- 
tributions have opposite signs and almost cancel. Furthermore, matrix element of individual 
operators contain three main components ("eights", "eyes", and "subtractions"), which again 



conspire to almost cancel each other (see Fig. [14]). Thus e'/e is extremely sensitive to each of 



these components, and in particular to their matching. 
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Figure 14: Three contributions to (Oq) matrix element: "eight" (boxes), "eye" (diamonds) and 
"subtraction" (crosses). These data represent bare operators for the dynamical ensemble. The 
fit is done for the sum total of all contributions. All errors were combined by jackknife. 



Consider fermion contractions with operators such asQ 

(PP)eu = (*75®&u)(u75®&t0 (33) 
(SS) W = (si <g> ld)(3l ® Id) (34) 
{PS) A 2U = {s-y 5 ®&d){dl®ld), (35) 

which are main parts of, correspondingly, "eight", "eye" and "subtraction" components of Oq 
and 0§ (see the Appendix). The finite renormalization coefficients for these operators have been 
computed in Ref. plf . The diagonal coefficients are very large, so the corresponding one-loop 



corrections are in the neighborhood of —100%. In addition, they strongly depend on which q* 
is used (refer to Table H). Thus perturbation theory fails in reliably matching the operators in 
Eqs. 



We apologize for slightly confusing notation: we use the same symbols here for operators as in the Appendix 
for types of contractions. 
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Table 2: (Oq) in arbitrary units with one-loop perturbative matching using two values of q* for 



Qi ensemble. 


For comparison. 


the results 


with no matchinj 


1 ("bare") 


are given. 


Quark mass 


0.01 


0.02 


0.03 


0.04 


0.05 


q* = 1/a 


0.1 ± 1.2 


-0.9 ±0.4 


-1.2 ±0.2 - 


-1.6 ±0.3 


-1.1 ±0.2 


q* = n/a 


-13.1 ± 1.8 


-9.0 ±0.5 


-7.1 ±0.3 - 


-6.3 ±0.5 


-4.6 ±0.5 


Bare 


-55.6 ±5.0 - 


-35.4 ±1.5 


-27.0 ±0.9 - 


22.3 ± 1.4 


-16.4 ± 1.5 



The finite coefficients for other (subdominant) operators, for example (PP)ef, {SS)eu 
and (SS)ef, are not known for formulation with gauge- invariant operators^. For illustration 
purposes, in Table |2] we have used coefficients for gauge non-invariant operators computed in 
Ref. [p4}| , but strictly speaking this is not justified. 



To summarize, perturbative matching does not work and some of the coefficients are even 
poorly known. A solution would be to use a non-perturbative matching procedure, such as 
described by Donini et al. [27]. We have not completed this procedure. Nevertheless, can we 



say anything about e'/e at this moment? 

5.3 Partially nonperturbative matching 

As a temporary solution, we have adopted a partially nonperturbative operator matching pro- 
cedure, which makes use of bilinear renormalization coefficients Zp and Z$. We compute the 
latter |2Bj following the non-perturbative method suggested by Martinelli et al. Namely 



we study the inverse of the ensemble-averaged quark propagator at large off-shell momenta in 
a fixed (Landau) gauge. An estimate of the renormalization of four-fermion operators can be 
obtained as follows. 

Consider renormalization of the pseudoscalar-pseudoscalar operator in Eq. (|3^). At one- 



loop level, the diagonal renormalization coefficient Cpp (involving diagrams shown in Fig. [L5|) 
is almost equal to twice the pseudoscalar bilinear correction Cp. This suggests that, at least at 
one-loop level, the renormalization of {PP) eu comes mostly from diagrams in which no gluon 
propagator crosses the vertical axis of the diagram (for example, diagram (a) in Fig. and 
very little from the rest of the diagrams (such as diagram (b) in Fig. |T3]). In other words, 
the renormalization of (PP)eu would be identical to the renormalization of product of two 
pseudoscalar bilinears, were it not for the diagrams of type (6), which give a subdominant 



5 Patel and Sharpe j24j have computed corrections for gauge-noninvariant operators. Operators in Eqs. ( |3£ 
( p5| ) have zero distances, so the corrections are the same for gauge invariant and non-invariant operators. Renor- 
malization of other operators (those having non-zero distances) differs from the gauge-noninvariant operators. 
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(a) 




(b) 



Figure 15: Example of one loop diagrams arising in renormalization of four-fermion operators: 
in type (a) no propagator crosses the axis, and type (b) includes the rest of the diagrams. 

contribution. Mathematically, 

(ppy E t = (PP) l EuZpp + ..., 

with 

Z PP = Z*(l + -^- 2 Cpp + 0{g 4 )) , (36) 

Zp = l + ^C P + Otf), (37) 
and dots indicate mixing with other operators (non-diagonal part). The factor Cpp = Cpp — 



2Cp contains diagrams of type (6) in Fig. [15] and is quite small. 

In order to proceed, it may be reasonable to assume that the same holds at all orders in 
perturbation theory, namely the diagrams of type (c) in Fig. |1^ give subdominant contribution 
compared to type (a) of the same Figure. This assumption should be verified separately by per- 
forming non-perturbative renormalization procedure for four-fermion operators. If this ansatz 
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is true, we can substitute the non-perturbative value of Z P into Eq. ( p6[ ) instead of using the 
perturbative expression from Eq. fl37|) . Thus a partially nonperturbative estimate of (PP)^ nt 
is obtained. This procedure is quite similar to the tadpole improvement idea: the bulk of di- 
agonal renormalization is calculated non-perturbatively, while the rest is reliably computed in 
perturbation theory. Analogously we obtain diagonal renormalization of operators (SS)ju and 
by using Z ss = Z%(\ + ^C ss + 0(g*)) and Z PS = Z S Z P (1 + ^C PS + 0(g*)). We 
note that Z P ^ Z s , even though they are equal in perturbation theory. We match operators 
at the scale q* = 1/a and use the continuum two-loop anomalous dimension to evolve to \i = 2 
GeV. 

Unfortunately, the above procedure does not solve completely the problem of operator renor- 
malization, since it deals only with diagonal renormalization of the zero-distance operators in 
Eqs. (|33] — [J^). Even though these operators are dominant in contributing to e'/e, other op- 
erators (such as (SS)eu and (PP)ef) can be important due to mixing with the dominant 
ones. The mixing coefficients for these operators are not known even in perturbation theory. 
For a reasonable estimate we use the coefficients obtained for gauge non-invariant operator 
mixing P^j . 

Secondly, since renormalization of operators (PP)eu, (SS)iu an d (PS)a2U is dramatic?], 
their influence on other operators through non-diagonal mixing is ambiguous at one-loop order, 
even if the mixing coefficients are known. The ambiguity is due to higher order diagrams (for 
example, those shown in Fig. In order to partially resum them we use operators (PP)eu, 
(SS)iu and (PS)a2U multiplied by factors Zp, Z\ and Z P Z S , correspondingly, whenever they 
appear in non-diagonal mixing with other operators []. This is equivalent to evaluating the 
diagrams of type (a) and (6) in Fig. [16] at all orders, but ignoring the rest of the diagrams (such 
as diagrams (c) and (d) in Fig |l^) at all orders higher than first. To estimate a possible error in 
this procedure we compare with a simpler one, whereby bare operators are used in non-diagonal 
corrections (i.e. we apply strictly one-loop renormalization). The difference in e'/e between 
these two approaches is of the same order or even less than the error due to uncertainties in 
determination of Z P and Z s (see Tables || and f|). 



6 For example, at m q = 0.01 and /i = 2 GeV for Q\ ensemble we obtain Zpp — 0.055 ± 0.007, Zps = 
0.088 ± 0.007 and Z ss = 0.142 ± 0.010. 

7 A completely analogous scheme was used for mixing of Oq with Oi through penguins when evaluating 
ReA . 
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Figure 16: Example of four kinds of diagrams with arbitrary number of loops arising in renor- 
malization of four-fermion operators: in (a) and (b) no propagator crosses the box or the axis; 
(c) and (d) exemplify the rest of the diagrams. The rectangular drawn in dotted line in (b) 
corresponds to operator structure PPeu- 
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m K 2 , GeV 2 

Figure 17: A rough estimate of e'/e for Q\ ensemble using the partially-nonperturbative proce- 
dure described in text. Three sets of points correspond to using experimental ReAo and Re/^ 
in Eq. ([TJD (crosses), using our ReA but experimental ui (diamonds), or using Rev4 and ReA 2 
obtained from our calculations (squares). All other details are the same as in Table |3|. The 
shown error is a combination of the statistical error, a matching error coming from uncertainties 
in the determination of Zp and Zg, and an uncertainty in non-diagonal mixing of subdominant 
operators. 



26 



6 Estimates of e'/e 



Within the procedure outlined in the previous section we have found that (Oq) has a different 
sign from the expected one. This translates into a negative or very slightly positive value of 
e'/e (Tables | and | and Fig. 0). 

If the assumptions about the subdominant diagrams made in the previous section are valid, 
our results would contradict the present experimental results, which favor a positive value of 
e' I e. They would also change the existing theoretical picture Jl2| due to the change of sign of 
(Oe). 

Finite volume and quenching effects were found small compared to noise. The main uncer- 
tainty in e'/e comes from operator matching, diagonal and non-diagonal. For diagonal matching 
the uncertainty comes from (1) errors in the determination of Zp and Z$ non-perturbatively 
and from (2) unknown degree of validity of our ansatz in Sec. |5.3| . For non-diagonal matching, 
the error is due to (3) unknown non-diagonal coefficients in mixing matrix and (4) ambiguity 
of accounting higher-order corrections. The error (1), as well as the statistical error, is quoted 
in Tables |3] and |j. The size of the error (4) can be judged by the spread in e'/e between two 
different approaches to higher-order corrections (strictly one-loop and partial resummation) , 
also presented in Tables |3| and [|. The error (3) is likely to be of the same order as the error 
(4). The error (2) is uncontrolled at this point, since it is difficult to rigorously check our 



assumption made in Sec. |5.3| . In Fig. 17 we combine the statistical error with errors (1) and 
(4) in quadrature. 

The uncertainty due to operator matching is common to any method to compute the relevant 
matrix elements on the lattice (at least, with staggered fermions). In addition, our method 
has an inherent uncertainty due to dropping the higher order chiral terms. Lattice spacing 
dependence of e'/e is unclear at this point, but it may be significant. 

We note that there are several ways to compute e'/e. One can use the experimental values 
of Re^o and Re^ in Eq. (|ToD, or one can use the values obtained on the lattice. One can 
also adopt an intermediate strategy of using the experimental amplitude ratio u> and computed 
KeA . When the higher-order chiral corrections are taken into account and the continuum limit 
is taken (so that uo = 22), these three methods should converge. At this point any of them can 
be used, and we compare them in Tables |] and £|. 

In view of the issues raised above, e'/e is an extremely fragile quantity. The rough estimates 
in Tables |3| and (| and Fig. [17| should be used with extreme caution. 
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7 Conclusion 



We have presented in detail the setup of our calculation of hadronic matrix elements of all 
operator in the basis defined in Eqs. (| — |ll]). We have obtained statistically significant data for 
all operators. Based on these data we make theoretical estimates of KeAo and Re^2 amplitudes 
as well as e'/e. 

Simulation results show that the enhancement of the AI = 1/2 transition is roughly consis- 
tent with the experimental findings. However, the uncertainty due to higher order chiral terms 
is very large. If these terms are calculated in the future, a more definite prediction for physical 
amplitudes can be obtained using our present data for matrix elements. Simulations should be 
repeated at a few more values of j3 in the future in order to take the continuum limit. 

Calculation of e'/e is further complicated by the failure of perturbation theory in oper- 
ator matching. We give our crude estimates, but in order to achieve real progress the full 
nonperturbative matching procedure should be performed. 

We appreciate L. Venkataraman's help in developing CRAY-T3E software. We thank 
the Ohio Supercomputing Center and National Energy Research Scientific Computing Cen- 
ter (NERSC) for the CRAY-T3E time. We thank Columbia University group for access to 
their dynamical configurations. 



Explicit expressions for fermion contractions 
.1 Quark operators 

We work in the two flavor traces formalism when calculating contractions with four-fermion 
operators: for each contraction separately the operators are rendered in the form (if necessary, 
by Fierz transformation) of two bilinears with the flavor flow in the form of a product of two 
flavor traces. To be more precise, for "eight" contractions the operators are rendered in the 
form (sTu)(uTd), while for the "eye" and "annihilation" contractions the appropriate form is 
(sTd) (qTq) . This is done in the continuum, before assigning the staggered fermion flavour. 
The operator transcription in flavor space for staggered fermions is now standard ||10|1 , and 



we give it here for completeness. The Goldstone bosons have spin-flavor structure 75 <g> £5. 
The flavor structure of the operators is defined by requiring non-vanishing of the flavor traces, 
and so it depends on the contraction type: the flavor structure is £ 5 in "eights" and two-point 
functions, 1 in "eyes" and "subtractions". In "annihilation" contractions the flavour structure 
is 1 for the bilinear in the quark loop trace and £5 for the one involved in the external trace. 
Either one or two color traces may be appropriate for a particular contraction with a given 
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operator (see the next Appendix section for details). In one trace contractions (type "F" for 
"fierzed" ) the color flow is exchanged between the bilinears, while in two trace contractions (type 
"U" for "unfierzed") the color flow is contained within each bilinear so that the contraction 
is the product of two color traces. In either contraction type, when the distance between 
staggered fermion fields being color-connected is non-zero, a gauge connector is inserted in the 
gauge-invariant fashion. The connector is computed as the average of products of gauge links 
along all shortest paths connecting the two sites. We also implement tadpole improvement by 
dividing each link in every gauge connector by u = (l/3Tr(Up)) 1 / 4 , where Up is the average 
plaquette value. 

.2 Sources and contractions 

We use local U(l) pseudofermion wall sources. Explicitly, we set up a field of U(l) phases 
£ a (x;i ) (|£a| = 1) f° r eac h color and each site at a given timeslice to, which are chosen at 
random and satisfy 

(C(x;t )^(y;t ))=^5 x , y . (38) 

(Boldface characters designate spatial parts of the 4- vector with the same name.) We proceed 
to explain how this setup works in the case of the two-point function calculation, with trivial 
generalization to "eight" and "annihilation" contractions. 

Consider the propagator from a wall at t = in a given background gauge configuration, 
computed by inverting the equation 

(p + m)f y xpiy) = 0)4 4 ,o. (39) 
This is equivalent to computing 

Xf>(y)=Y,Ux> Q ) G ^(y; x,o), (40) 

X 

where G(y; x) is the propagator from 4-point x to 4-point y. For staggered fermions description 
we label the fields by hypercube index h and the hypercube corner indices e {0, l} 4 instead 
of y. The two-point function is constructed as follows: 

TP ex X* a (h, A)U afS (h, A,A + A) X p(h, A + A) 0(A) (-1) A , (41) 

h,A 

where 4>{A) and A M are phases and distances appropriate for a given staggered fermion operator 
f] , U(h, A, A + A) is the appropriate gauge connector (see below), modulo 2 summation is 

8 For a given bilinear with spin-flavor structure Ts <8> IV, these are determined as follows: A M = |5 M — F^\ 2 
and 4>(A) = \Ty{t\ T s T a+a t\,), where and are spin and flavor vectors such that T s = 7^7^73 3 74 4 
and T F = 7 f 1 7 2 F2 7 3 F3 74 F4 , and T A = ^l^l^lt' ■ 
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implied for hypercube indices A, and h runs over all hypercubes in a given timeslice t where 
the operator is inserted. The factor (— 1) A takes into account that for staggered fermions 
G(x;y) = G^(y] x)(— l) x (— l) y - Equation ([IT]) corresponds to 

TP oc G a p(z,y) T G^y,x) (-1) 2 C a (z%(x) , (42) 

x,y,z 

where V is used for simplicity to show the appropriate operator structure. The summation over 
x and z over the entire spatial volume averages over the noise, so the last equation is equivalent 
to 

TP oc Y, to G{x, y) T G{y, x) (-1)*. (43) 

Therefore, using the pseudofermion wall source is equivalent to summation of contractions 
obtained with independent local delta-function sources. Note that the factor (— l) x and zero 
distance in the staggered fermions language are equivalent to spin-flavor structure 75 (g> £5. 
This means the source creates pseudoscalar mesons at rest, which includes Goldstone bosons. 
Strictly speaking, this source also creates mesons with spin-flavor structure 7574 ® £ 5 £ 4 , since 
it is defined only on one timeslice. However, as explained in the first footnote in Section 2.3 of 
Ref. jl(J, these states do not contribute. 



We have used one copy of pseudofermion sources per configuration. 

Analogously, we construct the pion sink at time T by using another set of U(l) random 
noise ((£*(x;T) <^(y;T)) = 5 a ^ <5 Xjy , |£| = 1). The propagator $ is computed as follows: 

(p + m)f y Mv) = £»(x; T)5 X4 , T . (44) 

Suppose Ax, A 2 G {0, l} 4 and <fii(A), 4>2(A) are distances and phases of the two staggered 
fermion bilinears making up a given four-fermion operator. The expression for the "eight" 
contraction (Fig. |TJa) with two color traces ( "U" type) is given by 

Eu« £ xl{h,A)U a p(h,A,A + A 1 ) X p{h,A + A 1 )MA)(-l) A 



h,A,B 



x$*(h, B)U pa (h, B,B + A 2 )<5> a (h, B + A 2 ) <f> 2 (B) (-1) B , (45) 



up to various normalization factors which cancel in the B ratio. In this expression A, B e 
{0, l} 4 run over 16 hypercube corners (modulo 2 summation is implied for these indices). The 
hypercube index h, as before, runs over the entire spatial volume of the timeslice t of the operator 
insertion. The gauge connector U(h, A, B) is the identity matrix when A = B, otherwise it is 
the average of products of gauge links in the given configuration along all shortest paths from A 
to B in a given hypercube h. The expression fl4"5|), as well as all other contractions, is computed 
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for each background gauge configuration and is subject to averaging over the configurations. 
(Whenever several contractions are combined in a single quantity, such as a B ratio, we use 
jackknife to estimate the statistical error). 

The expression for one color trace ( "F" type) contraction is similar: 

E F cx X* a (h,A)U a p(h,A,B + A 2 ) X „(h,A + A 1 )MA) (-1) A 

h,A,B 

x$* p (h, B)U pa (h, B,A + A 1 )^ p (h, B + A 2 ) fa(B) (-1) B , (46) 

For "eye" and "subtraction" diagrams (Fig. |l]b and |l|d) the source setup is a little more 
involved, since the kaon and pion are directly connected by a propagator. In order to construct 
a wall source we need to compute the product 

4>(v) = E G (y » *; ^ T ) ■ T ; o. °) 

X 

In order to avoid computing propagators from every point x at the timeslice T, we first compute 
propagator G(x, T;0,0), cut out the timeslice T and use it as the source for calculating the 
propagator to (y, £). This amounts to inverting equation 

(P + m)f y Mv) = Xa(x) 5 {X4 , T) (-1T , (47) 

where Xa( x ) is the propagator from the wall source at t = defined in Eq. (ftSQ. We use the 
following expression for evaluating the "subtraction" diagram: 

S oc XUK A)U aP (h, A,A + A)Mh, A + A) 0(A) (-1) A , (48) 

h,A 

Again, averaging over the noise leaves only local connections in both sources, so in the contin- 
uum language we get: 

s ex J2 tr G ( x » o; z > *) r G & *; y> T ) G (y- T ; x > °) • ( 4 9) 

x,y,z 

(In fact, we are mostly interested in subtracting the operator si <g> Id, so in Eq. (ffl) A = 
{0,0,0,0} and <p(A) = 1.) 

In order to efficiently compute fermion loops for "eye" and "annihilation" diagrams (Fig. |l]b 
and |l|c), we use U(l) noise copies Q % \ i = 1, . . . , N, at every point in space-time. We compute 
77 W by inverting (^) + m)r)( l > = Q l >. It is easy to convince oneself that the propagator from y 
to x equals 

G(x;y) = ( Vx C). (50) 
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In practice we average over N = 10 noise copies. This includes 2 or 4 copies of the lattice in 
time extension, so the real number of noise copies is 20 or 40, with another factor of 3 for color. 
The efficiency of this method is crucial for obtaining good statistical precision. 
The expression for "U" and "F" type "eye" diagrams are as follows: 



lu oc £ Xa(h,^W a p(h,A,A + A 1 )MKA + A 1 )<f> 1 {A)(-iy 
i 

1 



h,A,B 
N 

x T?E (^(h,B)U pa (h,B,B + A 2 ) V ^(h,B + A 2 )MB)(-l) B , (51) 

If oc £ x* a (h,A)U a(7 (h,A,B + A 2 )Mh,A + A 1 ) ( j )l (A)(-l) A 

h,A,B 
1 N 

x TtH Cf*{KB)U pP {h ) B,A + A 1 )4\h ) B + A 2 )^{B){-l) B . (52) 

i=i 

The computation of "annihilation" diagrams (Fig. |l|c) is similar to the two-point function, 
except the fermion loop is added and the derivative with respect to the quark mass difference 
ma — m s is inserted in turn in every strange quark propagator. Derivatives of the propagators 
are given by inverting equations 

(P + m) X ' = X, (53) 
(p + m)t/ w = r] (i) . (54) 

We have, therefore, four kinds of "annihilation" contractions, which should be combined in an 
appropriate way for each operator depending on the quark flavor structure (this is spelled out 
in the next Appendix section): 



Aiu « E Xa{h,A)U aP {h,A,A + A 1 ) X p(h,A + A 1 )MA){-l 
1 



h,A,B 
N 



x ^E C?*(h,B)U pa (h,B,B + A 2 )^\h,B + A 2 )4 >2 (B)(-l) B , (55) 



^=i 



A if «E Xa(h,A)U atr (h,A,B + A 2 ) X p(h,A + A 1 )MA) (~1) A 

h,A,B 

1 N 

x T?E (?*(h,B)U p ^h,B,A + A 1 ) V ^(h,B + A 2 )MB)hl) B , (56) 

JV i=l 

A 2U oc X* a (h, A)U aP (h, A,A + A x ) X p{K A + A 1 ) <j> x {A) (-1) A 

h,A,B 
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1 N 

x T?£ Cr( /i ' 5 ) C/ p-( /i ^' S + A ^ ) ( /i ^ + A 2)02(5)(-l) B ) (57) 
iv 1=1 

A 2F oc X* a (K A)U aa (h, A,B + A 2 ) X p(h, A + A,) <pi(A) (-1) A 

h,A,B 
1 N 

x 77 E (?*(h,B)U p ^h,B,A + A 1 ) V '^(h,B + A 2 )MB)(-l) B . (58) 

JV i=l 

Explicit expressions for matrix elements in terms of fermion 
contractions. 

Operators in Eqs. (Q — |TT|) can be decomposed into 1 = and 1 = 2 parts, which contribute, 
correspondingly, to AI = 1/2 and A J = 3/2 transitions. Here we give the expressions for these 
parts for completeness, since ReA , ReA 2 and e'/e are directly expressible in terms of their 
matrix elements. The 1 = parts are given as follows: 

Of = ^(^(l-75)^)(^(l-75)«)-^(^(l-75)«)(^ t (l-75)rf) 

+ l(s 7At (l- 75 )d)(d 7 ^(l- 75 )d) (59) 

Of = ^(^(l- 75 ) M )(^(l- 75 )d)-^(^(l- 75 )d)(nr(l- 75 ) M ) 



1 

+ 3 

)(0) 



(37^(1 - 75 )d)(dy(l-75)d) (60) 



0£> = (^(l- 75 )d) £ (?7 M (l-75)?) (61) 
Of = (^(1-75)^) E (^7 M (l-7 5 )g a ) (62) 
Of = (3^(1 - 7B )d) £ (g 7 ^(l+ 75 )g) (63) 

<J=lt, cf,s 

Of = (S a7M (l-7B)d/j) E (^7 M (l + 75)g«) (64) 

Of = i[(57M( 1 -T5)tO(^(l + 75)«)-(57 M (l-7B)«)(S7 M (l + 7B)tO 

- (s 7M (l- 75 )d)(l 7 ^l + 75 )s)] (65) 

Of = -[(S a7t (l - 75)^) (%7 M (1 + 75)Ma) - (S a 7^(l - 7b)W/9) (W/37 M (1 + 7s)d a ) 
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- (s a7 „(l-75)^)(s/37 M (l+75)sa)] (66) 

{ 9 ] = ~[(*7„(1 - l 5 )d)(uY(l - 7s)«) - - 7 5 )«)(^(1 - 75)d) 

- (s 7M (l- 75 )rf)(^(l-75)s)] (67) 

Oi? = i[(s7 M (l- 75 )«)(^(l-75)rf)-(^(l-75)rf)(«7 At (l-75)n) 

- (^(1-75)^(^(1-75)5)] (68) 

Expressions for 1 = 2 parts follows: 

Of } = Of = 2 -Of ] = 2 -0 { $ = ^(1 - 7 5 )«)(«7 M (1 - 7fe)d) 

+(3ry M (l - 75)d)(wy(l - 7b )«) - (s 7m (1 - 75 )d)(d7"(l - 7B )d)] (69) 

Of = ^[(^(l- 75 V)(«7 M (l + 75)rf) 

+(57 A4 (1 - 75 )rf)(ny t (l + 7 5 )w) - (^(1 - 75)rf)(3 7 M (l + 7 s)d)] (70) 

Of } = ~[(s a 7„(l- 75W)(%7 M (1 + 7 5 K) 

+ (s q7m (1 - 75)^)^7^(1 + 75K) - (s a7|4 (l - 75 )^)(^(l + 75 )d a )] (71) 

0f ] = of } = Of } = Of } = (72) 

(Whenever the color indices are not shown, they are contracted within each bilinear, i.e. there 
are two color traces.) 

As mentioned in Sec. |3.2| , in order to compute matrix elements of I = operators one 
needs to evaluate three types of diagrams: "eight" (Fig. |I|a), "eye" (Fig. [L|b) and "annihilation" 
(Fig. ^c). In the previous Appendix section we have given detailed expressions for compu- 
tation of these contractions, given the spin-flavor structure. Here we assign this structure to 
all contractions required for each operator, i.e. we express each matrix element in terms of 
contractions which were "built" in the previous section. 

Let us introduce some notation. Matrix element of the above operators have three compo- 
nents: 



(tt+tt-IO^ ) = (Ei + Ii-S ( 2 ma,) ) ™ K ™* (73) 



2 



where m is the common quark mass for s, d and u, and 

p 



Oi = -4- (74) 
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Here E{ and Ii stand for "eight and "eye" contractions of the (7r + \Oi\K + ) matrix element, 
A{ ~ (0\Oi\K°) I (md — m s ) is the "annihilation" diagram, S = (ir + \sd\K + ) is the "subtraction" 
diagram, and P = (0\s , ~f5d\K ) is the two-point function. We compute ai by averaging the ratio 
in the right-hand side of Eq. ([74]) over a suitable time range. 

Detailed expressions for Ei, Ii and Ai are given below in terms of the basic contractions 
on the lattice. We label basic contractions by two letters, each representing a bilinear. For 
example, PP stands for contraction of the operator with spin structure (75) (75), SS is for 
(1)(1), VV stands for (7^X7^), and AA is for (7^75) (7^75). The staggered flavor is determined 
by the type of contraction, as explained in the previous Appendix section. Basic contractions 
are also labeled by their subscript. The first letter indicates whether it is an "eight", "eye" 
or "annihilation" contraction, and the second is "U" for two, or "F" for one color trace. For 
example: PPeu stands for the "eight" contraction of the operator with spin-flavor structure 
(75 £$£5) (75 ® £5) with two color traces; VAaif stands for the "annihilation" contraction of the 
first type, in which the derivative is taken with respect to quark mass on the external leg (see 
the previous Appendix section), the spin-flavor structure is (7^ ® £5) (7^75 <S> 1), and one color 
trace is taken. What follows are the full expressions 9 . 

"Eight" parts: 



E? ] 


= 2 -{VV EF + AA EF )- l -{VV EU + AA EU ) 




(75) 




= \{VV EU + AA EU )-^{W EF + AA EF ) 




(76) 


eT 
eP 


= VV EF + AA EF 
= VV EU + AA EU 
= 2{PP EF — SS EF ) 
= 2{PP EU — SS E u) 




(77) 
(78) 
(79) 
(80) 


Ei 0) 


= SS EF - PP EF + l - (VV EU - AA EU ) 




(81) 




= SS EU -PP EU + ^(VV EF -AA EF ) 




(82) 


4 0) 


= = \iVV EF + AA EF - VV EU - AA EU ) 




(83) 


E[ 2) 


= E® = |^ 2) = |sg } = ^(VVeu + AA EU + VV EF + AA EF ) 


(84) 



9 Signs of operators O7 and 0% have been changed in order to be consistent with the sign convention of Buras 
et al. |12| . 
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4 2) 






(85) 


Ef 


= i(AA EU -VV EU ) + SS EF 


— PPef 


(86) 


Ef 


= ±(AA EF -VV EF ) + SS E u 


— PPeu 


(87) 



"Eye" parts: 

I? = VV W + AA W (88) 

/(°) = VVif + AAjp (89) 

/f = 3(VV W + AA W ) + 2(VV IF + AA IF ) (90) 

= 3(VV IF + AA IF ) + 2(VV IU + AA IU ) (91) 

4°) = Z(VV W -AA W )+A{PP IF -SS IF ) (92) 

4°) = ^VViF-AAj^+AiPPju-SSw) (93) 

4 0) = 2(PP IF -SS IF ) (94) 

4 0) = 2{PP IU -SS IU ) (95) 

4°) = VVip + AAj F (96) 

4? = V^ + A4 ro (97) 

"Annihilation" parts are obtained by inserting the derivative with respect to (m^ — m s ) into 
every propagator involving the strange quark: 

4°) = -(yA Alu + AV Alu ) (98) 

4 0) = -(VVW + AIW) (99) 
4 0) = -3(VA A1U + AV A1U )-(VA A2U + AV A2U ) 

-2(VA A1F + AV A1F ) - (VA A2F + AV A2F ) (100) 
4 0) = — 3(T^A^ijt + AV A \ F ) — (V A A2F + AV A2F ) 

-2{VA AW + AVk ll7 ) - (VA A2U + AV A2U ) (101) 

4°) = 3(VA A1U -AV A1U ) + (VA A2U -AV A2U ) + 2(PS A2F -SP A2F ) (102) 

4 0) = 3(VA A1F -AV A1F ) + (VA A2F -AV A2F ) + 2(PS A2U -SP A2U ) (103) 

4 0) = ^(^A2 C /-^V A2[/ ) + (P^ 2F - ( 5P A2i ,) (104) 

4 0) = \{VA A2F -AV A2F ) + (PS A2U -SP A2U ) (105) 
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(106) 



(107) 
(108) 



Of course, "eye" and "annihilation" contractions are not present in I = 2 operators. 
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TABLES 



Table 3: e'/e in units of 10~ 4 for Qi ensemble, computed in three ways: (Ml) Rev4 and KeA 2 
values entering the expression for e'/e are taken from experiment; (M2) u amplitude ratio 
is taken from experiment, while KeAo amplitude is from our simulation; (M3) both KeAo and 
KeA2 are from our simulation. Partially-nonperturbative matching have been used to obtain the 
results. In all perturbative corrections we have used one-loop coefficients computed for gauge- 
noninvariant operators. The first error is statistical (obtained by combining the individual 
errors in matrix elements by jackknife). The second error is the diagonal operator matching 
error due to uncertainty in the determination of Zp and Z$- In order to estimate the non- 
diagonal matching error we compare two renormalization procedures: using strictly one-loop 
non-diagonal corrections (denoted "(11.)"), and resumming part of higher-order corrections in 
non-diagonal mixing by using non-perturbative renormalization factors Zp and Z$ (as explained 
in Section |5T3| ). The latter method is denoted "(p.r.)" . Some other parameters used in obtaining 
these results are: ImA 4 = 1.5 -10~ 4 , m t = 170 GeV, m;, = 4.5 GeV, m c = 1.3 GeV, Q v + V > = 0.25, 
a^~°\2 GeV) = 0.195 (the latter is based on setting the lattice scale by p meson mass). 
Short distance coefficients were obtained by two-loop running using the anomalous dimension 
and threshold matrices computed by Buras et al. [12]. 



Quark mass 


0.01 


0.02 


0.03 


Ml (p.r.) 


-58.1 ±2.1 ±10.6 


-26.3 ±0.8 ±8.9 


-16.1 ±0.4 ±8.0 


Ml (11.) 


-52.3 ±2.2 ±10 


-22.0 ±0.8 ±8.3 


-12.2 ±0.5 ±6.9 


M2 (p.r.) 


-38.6 ±2.1 ±6.0 


-18.7 ±0.3 ±7.0 


-11.7±0.2±6.0 


M2 (11.) 


-45.4 ±3.5 ±8.6 


-18.8 ±0.4 ±7.0 


-10.3 ±0.3 ±6.0 


M3 (p.r.) 


-97 ± 14 ± 13 


-81 ±4 ±23 


-79±2±27 


M3 (11.) 


-142 ±28 ±29 


-88 ± 5 ± 35 


-75 ± 2 ± 39 


Quark mass 


0.04 


0.05 




Ml (p.r.) 


-8.0 ±0.9 ±7.2 


-4.4 ±0.9 ±7.2 




Ml (11.) 


-4.2 ± 1.1 ±6.5 


-1.2 ± 1.0 ±6.6 




M2 (p.r.) 


-6.1 ±0.5 ±5.3 


-3.1 ±0.5 ±4.9 




M2 (11.) 


-3.7 ±0.8 ±5.8 


-0.9 ±0.8 ±5.2 




M3 (p.r.) 


-81 ±5 ±38 


-74 ± 4 ± 38 




M3 (11.) 


-64 ± 4 ± 52 


-55 ±5 ±51 
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Table 4: e'/e results for Q3 ensemble (/? = 6.2). Everything else is the same as in Table |[ 



Quark mass 


0.010 


0.015 


Ml (p.r.) 


-109 ± 15 ±84 


-64 ± 7 ± 63 


Ml (11.) 


-95 ± 16 ± 74 


-52 ±8 ±56 


M2 (p.r.) 


-39 ± 9 ± 22 


-22 ±2 ± 18 


M2 (11.) 


-48 ± 18 ± 33 


-23 ±3 ±23 


M3 (p.r.) 


-30.0 ± 17 ±23.2 


-21.1 ±2.3 ±20.7 


M3 (11.) 


-48 ± 42 ± 39 


-24 ± 7 ± 28 


Quark mass 


0.020 


0.030 


Ml (p.r.) 


-36 ± 4 ± 60 


-19 ±2 ±38 


Ml (11.) 


-27±5±54 


-14 ±2 ±34 


M2 (p.r.) 


-13.7 ±0.6 ±20 


-9.5 ±0.6 ±17 


M2 (11.) 


-12.4 ± 1.1 ±23.3 


-8.1 ±0.9 ±18.5 


M3 (p.r.) 


-15.9 ± 1.8 ±26.5 


-10.8 ±1.1 ±21.6 


M3 (11.) 


-14.7 ± 1.4 ±31.1 


-9.6 ±0.8 ±18.4 
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